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The resilient modulus of different pavement materials is one 
of the most important parameters for the pavement design 
using the mechanistic-empirical (M-E) method. The resilient 
modulus is generally determined by a triaxial test, which is 
expensive and _ time-consuming and_ requires special 
laboratory facilities. This study aims to develop a model 
based on the Gaussian Process Regression (GPR) to predict 
the resilient modulus of stabilized base material with 
different additives under wetting-drying cycles. For this 
purpose, a laboratory dataset containing 704 records have 
been used. The input parameters were considered as_ the 
wetting-drying cycles, free lime to silica ratio, Alumina and 
iron oxide compounds in the additives, maximum dry density 
to optimum moisture content ratio, deviator stress, and 
confining stress. The results indicate high accuracy of the 
GPR method with a regression coefficient of 0.997 and 0.986 
respectively for train and test data and 0.995 for all datasets. 
Comparing the developed model based on the GPR method 
with the developed models in the literature based on the 
artificial neural network methods and the support vector 
machines shows higher accuracy of the GPR method. 
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1. Introduction 


The resilient modulus (M,) is measured to determine the material stiffness values within different 
stress levels and describe the stress-strain nonlinear behavior of soil under repeated loading. This 
parameter is one of the important key parameters for pavement analysis and design. In 1986, the 
AASHTO pavement design guide recommended using the M, to describe the stiffness of 
subgrade soil [1]. The M, has been used as one of the fundamental characteristics to describe 
materials in structural design of flexible pavement since the publication of the AASHTO 
guideline in 1986 [2]. Moreover, the resilient modulus is also used to describe soil characteristics 
and aggregate materials in the mechanistic-empirical pavement design methods (e.g MEPDG 
method) [3]. 


The resilient modulus of soil and the aggregate materials is defined as the applied deviator stress 
to recoverable axial strain ratio under dynamic loading as follows [4]: 


M = 24 (1) 


Where oy is the deviator stress and ¢, is the recoverable axial strain. In the uniaxial compression 
test without confined stress, og equals to the axial stress. 


The resilient modulus is directly determined in the laboratory by dynamic triaxial test, resonance 
column test, torsional shear test, and gyratory methods [5—7] or indirectly through correlation 
with the results of other standardized tests or backcalculation [8]. The most common method to 
determine the resilient modulus of soil and the aggregate materials in the laboratory is the 
dynamic triaxial test under the effect of different confining and deviator stresses. The dynamic 
triaxial test is time-consuming and expensive therefore, it is useful to present new methods and 
techniques to obtain an accurate estimation of the resilient modulus. Different statistical and 
computational intelligence models have been developed to predict the resilient modulus of fine- 
grained and coarse-grained soil so far, including the models developed based on the artificial 
neural network [9-13], support vector machines [14—16], adaptive neuro-fuzzy inference system 
[17] as well as hybrid models [18]. 


The stabilized road materials should be sufficiently durable under the effect of Wetting and 
Drying (W-D) as well as Freezing and Thawing (F-T) cycles [19]. According to the mechanistic- 
empirical pavement design guidelines, among other factors, the W-D and F-T cycles are 
important factors which cause damage and deterioration of base and subbase materials and likely 
lead to premature failure of the pavement [20]. Several researches have been carried out 
concerning the evaluation of the effect of W-D and F-T cycles on the mechanical and physical 
properties of the stabilized base layer [21—26]. In this regard, George and Davison (1963) carried 
out an experimental study to evaluate the durability of stabilized fine-grained soils under F-T 
cycles. The results of their study showed that applying ten F-T cycles is sufficiently destructive 
for the stabilized base layer [27]. Nonan and Hamfri (1990) evaluated the effect of the F-T cycles 
on the resilient modulus of the limestone aggregates stabilized with cement. Their results 
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indicated the importance of the F-T cycles over the W-D cycles at low deviator stresses. 
Moreover, they found out that the resilient modulus of the stabilized aggregates decreases 
significantly during the primary stages of the F-T cycles [21]. To evaluate the stabilized soils in 
Oklahoma, Miller et al. (2000) subjected stabilized samples to W-D cycles according to the 
ASTM D559 test method. In this study, the unconfined compressive strength (UCS) test was 
employed to evaluate durability. Results showed that the UCS increases by the increase in the 
W-D cycles [28]. Khoury and Zaman (2002) evaluated the effect of the W-D cycles on low 
quality aggregates stabilized with class C fly ash. Aggregates with low quality are referred to as 
the final acceptable limit to be used in the base layer in Oklahoma. In their study, the resilient 
modulus, the UCS test, and the elastic modulus were used to evaluate the effect of W-D cycles. 
The results showed that the resilient modulus values of the cured specimens in 28 days under 30 
W-D cycles are approximately 5% lower than the resilient modulus of similar specimens without 
exposure to the W-D cycles [29]. By evaluating the effect of F-T cycles on the flexural 
characteristics of the aggregates stabilized with 10% of class C fly ash, Khoury and Zaman 
(2000) found out that the resilient modulus and rupture modulus decrease as the F-T cycles 
increase. Based on the results, they noted that the effect of F-T cycles on these two mechanical 
properties is a function of curing time and the number of F-T cycles [30]. They also explored the 
effect of the W-D cycles on the resilient modulus of the aggregate materials stabilized with 
different stabilizers in 2007 [31]. The results showed that the changes in the resilient modulus 
under W-D cycles can be determined by the rate of the chemical reaction. It was also stated that 
changing the values of the resilient modulus can be expressed more properly by lime content, 
SFA (Si02+A1,03+Fe,203) content, the optimum moisture content, and the maximum dry density. 
These researchers also proposed a regression model to predict the effect of these parameters as 
well as W-D cycles on the resilient modulus changes [31]. Maloof et al (2012) used the support 
vector machine (SVM) to predict the resilient modulus of the aggregate stabilized materials 
subjected to W-D cycles. They showed that the support vector machine results in more accurate 
modeling in comparison with the Least Square (LS) method [20]. Ghanizadeh and Rahravan 
(2016) used the artificial neural network (ANN) to predict the resilient modulus of the stabilized 
base materials under W-D cycles and compared their results with Maloof et al (2012) study. 
Results of this study confirms higher accuracy of ANN in comparison with the SVM method [9]. 
Kaloop et al (2019) developed some models to predict the resilient modulus of the stabilized 
base materials by combining the particle swarm optimization algorithm with the artificial neural 
network (ANN-PSO) and the extreme learning machine (PSO-ELM) and concluded that the 
PSO-ELM method has higher accuracy compared to other methods [32]. 


This study aims to present a model based on the Gaussian process regression (GPR) to predict 
the resilient modulus of stabilized base materials subjected to W-D cycles. The developed model 
is compared to the models developed by other researchers. The importance of each of the 
parameters to predict the resilient modulus of the stabilized base materials is considered by the 
sensitivity analysis. Also, the effect of each of the input parameters on the resilient modulus is 
evaluated by parametric analysis. 
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2. Dataset 


In this research, the dataset was adopted from Maloof et al. (2012) study [20]. In this dataset, 
four types of aggregate including Meridian, Richard Sput, Sawyer, and Rhyolite were stabilized 
by different additives and tested under different W-D cycles to determine the resilient modulus. 
The Meridian is a limestone material with 97% calcium carbonate (CaCO3) and Richard Spur is 
a limestone material with 87% calcium carbonate. Also, the Sawyer aggregate is a type of 
sandstone having a SiO» amount of about 94% [20]. 


The stabilizer agents in Maloof et al. study included Cement Kiln Dust (CKD), Class C Fly ash 
(CFA), and Fluidized Bed Ash (FBA). The specimens were made near the optimum moisture 
content and maximum dry density and were cured within 28 days with an approximate 
temperature of 21°C and humidity percentage was about 90%. Then the specimens were exposed 
to 0, 8, 16, and 30 W-D cycles, and the resilient modulus of the stabilized materials were 
obtained to evaluate the performance of the stabilized materials subjected to W-D cycles. The 
resilient modulus test was carried out by applying a haversine dynamic load with loading time of 
0.1s and rest period of 0.9s [20]. Finally, a dataset containing 704 records, five input variables 
and one output variable was stablished [20]. 


Previous studies show that the resilient modulus of the stabilized aggregate base is a function of 
the Wetting-Drying Cycles (WDC), Calcareous/Siliceous Fly Ash Ratio (CSFAR), the dry 
Density to Moisture content Ratio (DMR), confining stress (03) and deviator stress (64) [31]. 
Therefore, in the current study these input variables are used to develop the GPR model. 


The statistical characteristics along with the frequency and the cumulative frequency for this 
dataset is shown in Figure 1. 


3. Modeling by the Gaussian process regression (GPR) 


3.1. Gaussian process regression (GPR) 


The Gaussian process regression is a probabilistic non-parametric learning method which is 
widely used for regression and classification problems [33]. This method has attracted the 
attention of many researchers in various scientific fields [34,35]. GPR is highly efficient for 
modeling the nonlinear data due to kernel functions. Moreover, the main merit of GPR is 
providing a reliable response for the input data [36]. 


Assume that in a training set of D= \(x,, y; )i = 1,....n}, X € R”™ is the input data points (design 
matrix) and y € R”is the desired output vector. The main assumption of GPR is that the output 
can be calculated as follows [37,38]: 


y=fte (2) 


Where ¢ ~ N (0, o, ) € R is the equal noise variance for all x; samples. 
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WDC: the number of Wetting-Drying Cycles 


CSAFR: Calcareous to Siliceous Fly Ash ratio (silica, Alumina, and iron oxide compounds in cementitious materials) 


DMR: Dry Density to Moisture content Ratio 
o3: Restrictive stress (kPa) 

od: Deviant stress (kPa) 

MR: resilient modulus (MPa) 


Fig. 1. The statistical characteristics, frequency and cumulative frequency of inputs and output variables. 


In the GPR method, n observation in y= {y,,...,y,} vector is considered as a single point 
instance of the Gaussian multivariate distribution. In addition, it can be assumed that this 
Gaussian distribution has a mean of zero. The covariance function k(x, x') determines the 
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relationship of one observation to another. The square of the exponential covariance function is 
commonly used to approximate functions using the GPR method, which is as follows [34,38]: 


1\2 
=H. 
kliaw)aa2% of bod + 026(x, x’) (3) 
Where the maximum allowable covariance is defined as Ors It is noteworthy that k(x, x’) is 
equal to the maximum allowable covariance only when x and x' are so close to each other, 
therefore, f (x) is approximately equal to f (a), Besides, | shows the length of the kernel 


function. Moreover, 0 (x, x') is Kronecker delta function which is defined as follows: 
6, =lifi=j and 6, =0 if i# j. 


Concerning the training dataset, the final objective of the learning process is to predict y« the 
output value for a new input pattern. To achieve this, it is essential to develop three covariance 
matrices as follows: 


According to the training data set, the ultimate goal of the learning process is to predict the 
output value of y« for a new input pattern. To achieve this goal, it is necessary to create three 
covariance matrices as follows: 


[ k(x,.%,) k(x,,2,) te k(x,,,,) ] 
k(x,,%,) k(x),x,) te k(x,x, ) 


| k(x, .%) ki) iad k(x.%,)| 


Ke=[eltv) Rte) kOe) (a 
K.. =k(x,,.x.) 


Regarding the assumption that the data are taken from a Gaussian multivariable distribution: 


pele el ° 


Since it has been proved that y, | y is developed from a Gaussian multivariate distribution with 
an average of K,K 'y and variance of K,,—K,K 'K: , the approximate average and variance of 


the predicted output are obtained as follows: 


Ely.)=K.K‘y 


6 
var(y.) = Kes —K.K 1K? e 
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After determining meta-parameters of kernel function, parameters of the model including k and 
oO, can be determined by Bayesian inference. After training, the GPR model can be used to 
predict unknown values by the known input values. 


3.2. Evaluation parameters of GPR model 


In this study, the RMSE (Root Mean Squared Error), MAE (Mean Absolute Error), MAPE 
(Mean Absolute Percentage Error), MAD (Mean Absolute Deviation), and R? (coefficient of 
determination) parameters have been used to evaluate the performance of GPR modeling. These 
parameters can be determined as follows: 


N 
RMSE = |-- 3°(m; — p; (7) 
N jt 
1 N 
MAE =— 3'|m; — p; 8 
wee Pi| (8) 
MAPE = + 9 —Pi)100 (9) 
N jz) Mj 
MAD = median ( |pi— median (p) ) (10) 
2 


Nd vn, - Ne iD N 
[nS ym -(©¥\m;) [vd pi Ly veih 


Where N is the number of data to evaluate the desired method, m; is the measured value for i® 
data point, and p; is the predicted value for i™ data point. 


4. Results and discussion 


For modeling by the GRP method, 80% and 20% of all dataset records have been considered as 
training and testing set, respectively. Training and testing sets were randomly selected from the 
total data points. For this purpose, the data were shuffled in the Excel program and then the first 
80% of the data points were selected as the training set and the second 20% of the data points 
were considered as a testing set. The linear basis function has been used for modeling by GPR 
method and model overfitting was controlled by considering 10% of the training set as the 
validation set. Exponential kernel function and linear basis function also were assumed for 
developing final models. These two functions was selected based on the try and error method and 
after evaluation of different functions. be In this study four GPR models including model “A” 
considering all 5 input variables, model “B” considering 4 input variables (by omitting oa), 
model “C” considering 4 input variables (by omitting 03), and model “D” considering 3 input 
variables (by omitting og and 03) were developed. Performance of these four GPR models is 
demonstrated in Figure 2 to 5. 
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As can be seen, sequence of four models developed in terms of accuracy (R” and RMSE) is as 
A>C>B>D. In fact, model “A” with 5 input variables and values of R’ and RMSE equal to 0.998 
and 5858.336 for the training set and 0.985 and 6427.365 for the testing set has the best 
performance and model “D” with 3 inputs and values of R? and RMSE equal to 0.922 and 
7631.497 for the training set and 0.912 and 1373.576 for the testing set has the lowest 
performance among the four developed models. 
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Fig. 2. The predicted versus measured resilient moduli based on the train and test datasets for model “A”. 
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Fig. 3. The predicted versus measured resilient moduli based on the train and test datasets for model “B”. 


Also omitting og variable from input parameters in model “B”, lead to more reduction in model 
accuracy compared to omitting the 03 variable in model “C” which confirms that the resilient 
modulus of the stabilized base material is more dependent on og parameter compared to the 63 


88 


parameter. The approximate equality of the coefficient of determination (R’) for the training and 
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testing set in all cases indicates the high generalizability of the developed models. 
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Fig. 4. The predicted versus measured resilient moduli based on the train and test datasets for model “C”. 
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Fig. 5. The predicted versus measured resilient moduli based on the train and test datasets for model “D”. 


5. Comparing of performance with other methods 


As mentioned in section 2, the dataset used in the present study was adopted from the results of 
Maloof et al. (2012) [20]. In their research, in addition to conducting laboratory studies, they 
presented a model for predicting the resilient modulus using the support vector machine (SVM) 
method. Comparison of the results obtained from the present study and Maloof et al. (2009) 
study is provided in Table 1. It is evident that the GPR method is more accurate than the support 
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vector machine method; for example, in case of Model “A”, the value of R’ for the whole dataset 
using the GPR method is 0.995, and for the SVM method is 0.6400. Table 1 also shows the 
results of modeling accuracy based on the PSO-ELM method based on the Kaloop et al. (2019) 
study [32]. As can be seen, depending on the number of input parameters, the accuracy of the 
GPR method can be more or less than the PSO-ELM method. However, if all input parameters 
are considered, the accuracy of the GPR model is much higher than the PSO-ELM model, so that 
the RMSE for the GPR model is less than half of RMSE for the PSO-ELM model. 


Table 1 
Comparison of the GRP method and previously developed models. 
GPR (This study) SVM [20] PSO-ELM [32] 
Model Inputs 
R RMSE R RMSE R? RMSE 
A WDC, CSAFR, DMR, 63, 64 0.995 116.894 0.640 1134592 0.981 253.439 
B WDC, CSAFR, DMR, 63 0.931 491.145 0.875 659.750 0.948 415.554 
C WDC, CSAFR, DMR, oa 0.968 329.743 0.959 371.309 0.973 304.451 
D WDC, CSAFR, DMR 0.920 525.932 0.900 593.540 0.921 521.080 


6. Sensitivity analysis 


In this study, the Cosine Amplitude Method (CAM) is employed to determine the degree of 
importance of each input parameter to predict the resilient modulus using the degree of 
sensitivity index. The degree of sensitivity index can be calculated by the following equation 
[39,40]: 


ya 
XV; 
2 j=l yd 
a n n 
2 2 
pS 129 Da jo) 
Da IS 


Where x; shows the i” independent variable for the j“ data point and y; shows the dependent 
variable for the Sie data point (respect to x,). To estimate the relationship between input and 
output variables, the value of R; must be close to 1, while R; with a value of zero virtually 
eliminates the possibility of extracting a correlation. Figure 6 shows the degree of importance of 
input variables based on the results from the measured and predicted values of the resilient 
moduli. As it can be seen, the importance of different parameters can be displayed as DMR > 

CSAFR > og > 63>WDC. In other words, the DMR is the most important parameter and the 
WDC is the least important parameter for predicting the resilient modulus of the stabilized base 
materials under W-D cycles. Also the difference of R; between the predicted and measured values 
for the resilient modulus for WDC, CSAFR, DMR, o3, and og parameters is 0.31%, 0.18%, 
0.19%, 0.22%, and 0.22%, respectively which shows the high accuracy of the GPR method to 
predict the resilient modulus. 


(12) 
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Fig. 6. The importance of each of the variables based on the CAM method. 
7. Parametric analysis 


Time and cost limitations as well as the limited access to appropriate equipment are generally the 
main obstacles in laboratory studies. In most cases, investigating the effect of each variable on 
the test results over a wide range of cases requires fabricating several specimens which are time- 
consuming and expensive. One of the advantages of modeling is to employ the developed 
models for parametric studies and evaluating the effect of each input parameter on the model 
output. As mentioned before, in this study the input parameters are WDC, DMR, CSFAR, 
confining stress (63), and deviator stress (64). 


In this study, by using the optimal GPR model, the effect of interaction between WDC and 
DMR, the interaction between DMR and CSFAR, the interaction between WDC and CSFAR, 
and the interaction between o3 and og on the resilient modulus of the stabilized base materials 
have been evaluated. For this purpose, the desired parameter was changed between its minimum 
value and its maximum value, and other parameters were considered equal to the mean values, 
and then the resilient modulus was determined according to the change of the desired parameter 
based on the optimal GPR Model. 
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Fig. 7. The parametric analysis of the effect of the input variables on the resilient modulus. 


As evident, the resilient modulus dramatically increases as the DMR increases. While the 
increase in W-D cycles and CSFAR slightly decrease and increase the resilient modulus, 
respectively. The increase in the confining and deviator stress leads to the increase in the resilient 
modulus and it can be seen that changes in the deviant stress have more effect on the resilient 
modulus than the restrictive stress. The interaction of confining and deviator stress shows that 
increasing these two values increases the resilient modulus and it can be seen that the changes of 
deviator stress have a greater effect on the resilient modulus than the confining stress. 


8. Conclusion 


In this study, the GPR modeling method was used to predict the resilient modulus of the 
stabilized base materials. The results of this research can be summarized as follows: 


1- With respect to the values of R’ and RMSE, model “A” with five inputs and values of R? 
and RMSE respectively equal to 0.998 and 85.574 for the training set and 0.985 and 
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242.624 for the testing set has the best performance and model “D” with three inputs and 
values of R* and RMSE respectively equal to 0.922 and 513.173 for training set and 0.912 
and 576.968 for testing set has the lowest performance. 


2- It was found that omitting og variables in model “B” lead to more reduction of model 


accuracy compared to omitting o3 in model “C” which shows a greater dependence of the 
resilient modulus of the stabilized base materials to the deviator stress compared to the 
confining stress. 


3- The results of sensitivity analysis show that the degree of importance of different input 


parameters on the resilient modulus is as DMR> CSAFR> og> 03> WDC. 


4- The difference of degree of sensitivity (R;) between the predicted and measured values of 


the resilient modulus for WDC, CSAFR, DMR, o3, and og parameters are 0.31%, 0.18%, 
0.19%, 0.22%, and 0.22% respectively. 


5- Results of the parametric analysis shows that the increase in DMR leads to an increase in 


the resilient modulus and the increase in the W-D cycles and CSFAR slightly decreases 
and increases the resilient modulus, respectively. It can also be seen that increasing the 
deviator stress and the confining stress will increase the resilient modulus of the stabilized 
base materials. 
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